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In the present work, we study dark solitons in dynamical lattices with the saturable nonlinearity 
and compare them with those in lattices with the cubic nonlinearity. This comparison has become 
especially relevant in light of recent experimental developments in the former context. The stability 
properties of the fundamental waves, for both on-site and inter-site modes, are examined analytically 
| and corroborated by numerical results. Furthermore, their dynamical evolution when they are found 

, to be unstable is obtained through appropriately crafted numerical experiments. 
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INTRODUCTION. 



In the past few years, the investigation of dynamical nonlinear lattices, where the "evolution variable" is continuum, 
while the "spatial variables" are discrete, has seen a tremendous growth. One of the main directions that has triggered 
both experimental and theoretical efforts in this context has been the development of optically induced lattices in 
photorefractive media, such as Strontium Barium Niobate (SBN) pj, and its experimental realization @, 0|- As a 
i— i. result, a remarkable set of nonlinear waves has been predicted and experimentally observed; these include (but are 
not limited to) discrete dipole 0, quadrupole 0, necklace @ and other multi-pulse/multi-pole localized patterns Q, 
impurity modes Q , discrete vortices Q , and solitons in radial lattices ^(j ■ A recent review of this direction can be 
■ found in |TT| . 

Another direction, also relevant to nonlinear optics, which has been growing in parallel, concerns the intriguing 
interplay of nonlinearity and discrete diffraction emerging in fabricated AlGaAs waveguide arrays |l2j| . A variety of 
phenomena, such as discrete diffraction, Peierls barriers, diffraction management and the formation of gap solitons 
I . |14j among others, have been experimentally traced in the latter setting. These phenomena, in turn, stimulated a 
' tremendous increase in the number of the theoretical studies concerned with such effectively discrete media (see, e.g., 
(T) , the reviews p^lll6| ). 

Additionally, in the past decade, there has been yet another context where the consideration of discrete lattices and 
nonlinear waves thereof are relevant; this is the soft-condensed matter physics of Bose-Einstein Condensates (BECs) 
|17|. Droplets of such condensates may be trapped in a periodic optical potential commonly known as "optical 
lattice" (see, e.g., and references therein). The physics of BECs confined in optical lattices has also experienced 
a tremendous growth over the past few years, leading to many important theoretical and experimental results; these 
include the prediction and manifestation of modulational instabilities [l^, the observation of gap solitons [2(| and 
Landau-Zener tunneling |2l| among many other features; reviews of the theoretical and experimental findings in this 
area have also recently appeared in [22I Elf . 

Many of the above studies of coherent structures (especially, as regards one-dimensional ones) have been conducted 
in self- focusing media (i.e., media with "attractive" nonlinearities), where bright solitons or multi-soliton versions 
thereof emerge. Dark solitons in lattices with the so-called self-defocusing (or "repulsive" , in the language of BECs) 
nonlinearity have, in some sense, been far less popular with only a few examples of relevant experimental studies. 
In particular, the first observation (to the best of our knowledge) of dark solitons arose in |2J| in the anomalous 
diffraction regime of AlGaAs waveguide arrays, which feature the Kerr-type cubic nonlinearity. On the other hand, 
very recently, dark solitons were realized in defocusing lithium niobate waveguide arrays, which exhibit a different 
type of nonlinearity, namely a saturable, defocusing one due to the photovoltaic effect |25| . 

The scope of the present work is to compare the above settings, namely that of the cubic nonlinearity with that of 
the saturable one, regarding the properties of dark solitons. In particular, our aim here is to examine both of these 
models and make inferences on the properties of the dark solitons on the basis of existence and stability analysis. 
We aim to mathematically justify the stability of the pertinent dark soliton modes (such as those centered on a site 
and those centered between two lattice sites) and to numerically corroborate our analytical findings. Furthermore, we 
intend to examine more quantitatively statements, such as those of [25J about improved stability of inter-site modes, 
and to examine whether phenomena similar to those emerging for bright solitons (such as the stability alternation 
between on-site and inter-site solutions |26() can occur for dark solitons. Our work also builds on earlier theoretical 
work on such structures which considered them principally in the cubic nonlinearity setting such as e.g., p^ . l28ll2^| . 
We expand the considerations of these earlier works by considering in more detail the stability problem and obtaining 
information about the key eigenvalues (and doing that systematically, starting from the anti-continuum limit, where 
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the lattice sites are uncoupled), as well as applying these ideas to a different setting, such as that of the saturable 
nonlincarity and readily comparing the results of the two cases. 

The paper is organized as follows. In the next section, we will present our analytical considerations both for the 
saturable and for the cubic case. Subsequently, we will focus on numerical results and examine how they compare with 
the analytical predictions. In the fourth section, we will summarize our findings, present our conclusions, and indicate 
some interesting directions for future studies. Finally, the appendix contains an alternative, perturbation-theory based 
description of our analytical results that we will illustrate to be equivalent (for small couplings) to the analysis of the 
earlier sections (and is included for reasons of completeness of the mathematical part of our exposition) . 



II. MODELS AND THEORETICAL SETUP 



A. Cubic Nonlinearity 



We will start the discussion of our analytical results by considering the somewhat simpler (at the level of math- 
ematical details) case of the lattice exhibiting the cubic nonlinearity. The relevant model is the discrete nonlinear 
Schrodinger (NLS) equation of the following dimensionlcss form: 

iu„ = -eA 2 u n + Pi\u n \ 2 u n , (1) 

where u n denotes the complex electric field envelope in the case of waveguide arrays, or the wavefunction of the 
BEC droplets in the wells of a deep optical lattice. The parameter e is the coupling constant that characterizes 
the tunneling coefficient, while the parameter f3i > (which can be scaled out but will be kept for the purposes 
of completeness/comparison with the saturable model) sets the strength of the nonlinearity. The overdot denotes 
differentiation with respect to the evolution variable t (which is the propagation distance in optics, and time in 
BECs), and, finally, A 2 u n — (u n +i + u n -i — 2u n ) stands for the discrete Laplacian. Notice that we will be, in 
particular, concerned with steady state solutions of the form 

u n = v n exp(— iAji), (2) 

where Ai denotes the propagation constant (or the chemical potential in the BEC context). This leads to the steady 
state equation: 

f3i\v n \ 2 v n - Axv n - eA 2 v n = 0. (3) 

We will also be concerned with the linear stability of the ensuing solutions which can be studied using the ansatz 

u„ = cxp(-iAit) [v n + $ (exp(Xt)p n + exp(X*t)q n )] , (4) 

where 5 is a formal linearization parameter and the resulting linear problem (to 0(<5)) for the eigenvalue A and the 
eigenvector (p„, q n ) T determines the stability of the configuration as follows: a configuration will be (neutrally) stable 
for this Hamiltonian system if VA, the real part A r of the eigenvalue (A = A r + zAj) is such that A r = 0. This is due to 
the fact that the Hamiltonian nature of the system enforces that if A is an eigenvalue, so are A*, —A and — A*. Note 
that in the above expressions * denotes complex conjugate, while T denotes transpose. The linear stability problem 
can be written in compact notation in the form: 

[p n \_[2p 1 \v n \ 2 -A 1 -eA 2 0iv* 
\Qn)~\ -AW* A 1 -2p 1 \v n \* + eA 2 




1. Anti- Continuum Limit 



The anti-continuum (AC) limit of the equation is characterized by e = 0. Then, the steady states are straightfor- 
wardly characterized by: 

/Al ,0jexp^„), (6) 

where 9 n is an arbitrary phase. We will restrain ourselves here (without loss of generality, see e.g., [3(j) to real 
solutions with 9 n S {0, 7r}. 
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We will be concerned with two types of dark soliton solutions, the well-known on-site and inter-site one (see e.g. 
113, The former one, at the AC-limit, is of the form: 

V n <-1 

v„= 

V„>1 

while the latter type of dark soliton reads: 

v n <o 



v n >i 

One can subsequently examine the linear stability of these prototypical configurations, as a starting point for the 
finite e case. In the AC limit, the linear stability simplifies greatly due to the fact that the sites become uncoupled. 
One can then easily see that for all non-zero sites, the relevant stability matrix of the AC-limit is identical and has 
the form: 

Aa ( \ \ ) (12) 

This yields a pair of zero eigenvalues for each of these non-zero sites. Hence, in an inter-site configuration at the 
AC-limit, the linearization would only result in such zero eigenvalues. 

The only difference of an on-site configuration lies in the existence of the central vq = site. The latter produces 
a 2 x 2 stability matrix of the form 

Ax ( - 1 J ) , (13) 

and, therefore, an eigenvalue pair A = ±iAi. 




2. Finite Coupling Case 

Let us first consider the solution profile. It is clear that for finite e the solutions will be deformed from their AC-limit 
profile of Eqs. Q~ CH - To address this deformation, we can expand the steady state solution into a power series: 

v n = v^+ev^+0(e 2 ). (14) 
The leading order correction can be straightforwardly computed by the substitution of the expansion into Eq. , as 

for all excited sites. For the zeroth site of the on-site configuration, the symmetry of the profile yields a zero correction 
(to all relevant orders). It is easy to see that the correction of Eq. I|15(l only contributes to leading order to the sites 
with n G {1, —1} for the on-site and to those with n € {0, 1} for the inter-site configuration. These corrections amount 
to: 

vi ] = i (16) 

- (17) 
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for the on-site and to: 



' (18) 



v£ ] = -v[ 1] (19) 

for the inter-site. 

The next step is to consider the stability problem. Clearly, the latter will have eigenvalue contributions from 
two sources. One will be the continuous spectrum that will emerge from the background (that corresponded to 
zero eigenvalues at the AC limit) and the other will be from the point spectrum in the vicinity of the center of the 
dark soliton configuration. We commence from the easier calculation of the former type. The continuous spectrum 
corresponds to plane wave eigenfunctions of the type {p n ,Qn} ~ exp(ifcn). These, in turn, result into a matrix 
eigenvalue problem, where the matrix is of the form: 

where si = Ai + 4esin (k/2). This leads to a continuous eigenvalue spectrum described by the dispersion relation 



u = iX = ±V s\ - Af, (21) 



which is associated with the eigenvalue band iX € [— y/l6e 2 + 8eAi, %/16e a + 8eAi]. 

We now turn to the point spectrum eigenvalues, stemming from the central part of the excitation. For the inter-site 
configuration, in order to address this issue, we consider to a leading order approximation only the two sites (n = 
and n=l) participating in the dark soliton (as they are the only ones modified to leading order in perturbation 
theory). Using the perturbative expansion of Eqs. (|10fl - Hll[) within the relevant part of the eigenvalue problem, we 
obtain the 4x4 matrix 



/ Ai - 2e -e Ai - 2e 

-e Ai - 2e Ax - 2e 

-Ai + 2e -Ai + 2e e 

\ -Ai + 2e e -Ai + 2e 



(22) 



which, importantly, leads to a pair of real eigenvalues 

A = ± v / 2eAi - 5e 2 , (23) 

which render the configuration immediately unstable, as soon as e becomes non-zero. This prediction will be compared 
with the numerical results in the following section. 

Finally, for the onsite configuration, one can use a similar argument entailing the three central sites of the solitary 
structure and constructing a 6 x 6 matrix whose eigenvalues can, in principle, be computed. However, the resulting 
expressions are too cumbersome to be useful, hence, here we offer two alternative arguments that, as we will show, 
describe very accurately the behavior of the relevant point spectrum eigenvalue originally at A = ±iAi, associated 
with the on-site soliton. 

The first approach that we offer is a rigorous one and is based on the so-called Gerschgorin's theorem [3l| . Let us 
consider matrices A — [aij] of order N and define the radii r; = j^i \o-ij\ and denote the circles in the complex 

spectral plane Zi — {z £ C : \z — au\ < ri}. Then, Gerschgorin's theorem states that the eigenvalues of the matrix 
belong to these circles and, in fact, its refined version states that if m of these circles form a connected set, S, disjoint 
from the remaining N — m circles, then exactly m eigenvalues are contained in S. Our problem is an excellent testbed 
for the application of Gerschgorin's theorem because the sole eigenvalue discussed above is at ±iAi in the AC-limit (of 
zero radius for the Gerschgorin circles), while all others are located at the origin. Hence, for small e <C Ai, the (single) 
relevant point spectrum eigenvalue remains in the corresponding Gerschgorin circle which can be easily computed. In 
fact, considering that only the diagonal and super- and sub-diagonal elements are present for a site with vq = 0, the 
Gerschgorin estimate for the relevant eigenvalue emerges immediately as: 

\iX ± (Ai - 2e) | < 2e (24) 

which, in turn, necessitates that the relevant eigenvalue lies between Ai — 4e < iX < Ax. This is a rigorous result 
based on the above theorem, which provides a definitive (linear) bound on the growth of the relevant eigenvalue. 
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On the other hand, there is also an even more successful approach (which, however, is to a certain degree an 
approximation), according to which, one considers the eigenvalue equations for the eigenvectors p n and g*. Considering 
the anti-symmetry of the central site one may use the approximation of p n +i + Pn-i ~ (in lines similar to [2?1 l28l|. 
which, however, postulated the stronger condition p n +i ~ p n — l ~ 0)- In the latter setting, the eigenvalue equations 
decouple and provide a relevant estimate for this eigenvalue as 

i\ = ± (At - 2c) . (25) 

This prediction will also be compared with our numerical results in the next section. It is also worth noting here that 
this eigenvalue, moving towards the spectral plane origin, will collide with the band edge of the continuous spectrum 
when the predictions of Eq. I|25|l and of the band edge of the continuous spectrum coincide, which occurs for 

O /o q 

e cr = Ai = 0.07735Ai. (26) 



This collision, per the opposite Krein signature 32, 33] of the relevant eigenvalues, will lead to a Hamiltonian Hopf 
bifurcation [34j and a quartet of eigenvalues, similarly to what was found in [27j (for relevant details, see below). 

B. Saturable Nonlinearity 

The second model that we will consider is a discrete dynamical lattice exhibiting a saturable defocusing nonlinearity 
|26l, l35| , which has recently been used in a variety of studies and pertains to the dark soliton experiments reported in 
[25j . In this case, the relevant discrete NLS equation is of the form: 

iu» = -eA 2 u„ - - f 2 u n , (27) 

1 + \UnV 

where the nonlinearity parameter j3 2 > 0. Similarly to the previous case, we use the standing wave ansatz of the form 
U n — exp(iA 2 f)i;„, to obtain the steady state equation 

/3 2 

A 2 f n - — — ; — r~v n - eA 2 v n = 0. (28) 

1 + Kr 

Using a similar ansatz to that of Eq. (0J, we then obtain the following stability equations: 

i\ f Pn I - I 2 l+l^l 2 + (1+KI 2 ) 2 6 2 (l+|f„| 2 ) 2 | / Pn \ (9Q \ 

-tS&w -*. + TrftF-n9fefi.+«A,JUJ- (29) 

1. Anti- Continuum Limit 

As before, we commence our analysis of the model at the AC-limit of e = 0, where the solutions can be found 
explicitly in the form 



^--l,olexp(i0 n ), (30) 



A 



2 



where we will assume, without loss of generality, that A 2 < /3 2 . As in the case of the cubic nonlinearity, the on-site 
dark soliton of the AC-limit has the form 



v n <-i = \Jy 2 - 1 (31) 

v n=0 = (32) 
«„>! = -^-1. (33) 
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while the inter-site one reads 



r„ „ = (34) 



«»>i = -y^' 1 (35) 

Now, we turn to the stability problem in the AC-limit. Once again, there are two types of blocks in the stability 
matrix, one associated with sites belonging to the background, and another related to the central zero site of the 
on-site configuration. The former 2x2 block is of the form 

leading once again to zero eigenvalues for the excited sites. On the other hand, for the site with vq — 0, the related 
block is: 

(A 2 -ft)(j \ (37) 
leading to an eigenvalue pair A = ±i(A2 — P2), which is again isolated from the rest of the (zero) eigenvalues. 



2. Finite Coupling Case 

For the finite coupling case, we once again expand the solution in a power series in e as in Eq. I|14|) and obtain the 
leading order correction corresponding to the excited sites as: 



(0) 



= (38) 

l+(i>i 0> ) 2 

This results into the leading order corrections for the on-site mode: 



while for the inter-site mode: 



(39) 



(40) 



(41) 



(42) 



Having examined the leading order perturbation to the solution, we now focus on the stability for the finite 
coupling case. The continuous spectrum can be found, as before, by assigning {p n , q n } ~ exp(ifcn); using the notation 
s% = (A|//32)(/32/A2 — 1), as well as $2 — S2 + 4esin 2 (fc/2), the relevant matrix becomes 



S2 S2 

-82 -h 



which yields the continuous spectrum 



(43) 



iX = ±J$l - si (44) 



and the phonon band iX € [— \/l6e' 2 + 8es 2 , \/16e' 2 + 8es2]. 

Finally, we consider the point spectrum associated with the inter-site and the on-site configurations. For the former, 
once again we consider the 4x4 matrix block associated with the two central sites. However, given the additional 
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complexity of the saturable case, in order to obtain an analytically tractable leading order expression, we make a 
further simplification here. Namely, in all the elements of the matrix that have 0(1) terms with 0(e) corrections, we 
only consider the 0(1) terms. This leads to the following matrix 



/ 


S2 


— e 


S2 







— 6 


S2 





S2 


V 


-32 





-82 


e 





S2 


e 


-82 



(45) 



which, importantly, leads to a pair of real eigenvalues 

A = ±y/2es 2 - e 2 , (46) 

that will be compared with our numerical findings. 

On the other hand, for the onsite case, similar considerations to those of the cubic case simplify the results consid- 
erably (despite the additional complications of the saturable case linearization). In particular, we can still apply the 
Gerschgorin argument to obtain a rigorous bound on the eigenvalue associated with the zero site as 

|*A ± (A a - jSa - 2e) | < 2e. (47) 

This results in the inequality 02 — A2 — 4e < iX < 02 ~ A2. 

However, again a very accurate description (as will be evinced by our numerical computations below) can be 
provided by the anti-symmetric approximation of p n +i + Pn-i ~ 0, which leads to the eigenvalue pair: 

iX = ± (0 2 - A 2 - 2e) . (48) 

Lastly, the collision of this pair with the continuous spectrum band derived above is what will produce the instability, 
hence equating the above expression with that of Eq. (|44|l . we obtain the critical coupling of 



y*2 + (/3 2 - A 2 )2 + (/3 2 - A 2 )s 2 - 2s 2 - (02 - A 2 ) 



(49) 



Before we move to our numerical results and their comparison with our analytical findings let us point out here 
that the main stability results obtained above will also be confirmed in the Appendix using a different mathematical 
technique. In particular, for completeness in the Appendix we consider a perturbative formulation of the stability 
problem and identify the key eigenvalues through their formal series expansion in powers of e (similarly to what was 
done above for the corrections to the solution itself). In this alternative way, we will show how to obtain the leading 
order analogs of Eqs. 1(23(1 [for cubic, inter-site], 1(25(1 [for cubic on-site], 146(1 [for saturable, inter-site] and 148(1 [for 
saturable, on-site]. 



III. NUMERICAL RESULTS AND COMPARISON 



In our numerical results we will consider the solutions to the two models as a function of the coupling constant e, 
for the following selection of parameters (/?i,Ai) = (1,1) for the cubic case and (/?2,A 2 ) = (1,0.5) for the saturable 
one. Recall that 0\ and /3 2 , the parameters characterizing the strength of the nonlinear terms, can always be scaled 
out with a rescaling of e and a re-parametrization of time. This indicates that they can be chosen to be equal to 
unity; this way, the two models also possess the same type of nonlinearity for small amplitude solutions (as can be 
readily observed by Taylor expanding the saturable model for |m| 2 <§; 1). In the cubic nonlinear model, Ai can also 
be absorbed in a rescaling of e and of the solution amplitude, hence we chose a typical value pertaining to a unity 
background (but our calculations above and some discussion below illustrate how the results depend on Ai). Given 
the above three choices, A 2 is then chosen in a way such that the two models have the same background intensity. 
Generally, the latter requirement imposes the condition 



Pi a 2 



1, (50) 



to which we return below. 

Firstly, in Fig. ^we show a typical comparison of the numerically obtained inter-site, as well as on-site exact (up to 
a prescribed numerical tolerance) solutions for a typical value of e (namely, e = 0.5). These results have been obtained 
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FIG. 1: (Color online) Profiles of the dark soliton solutions of the two models for unity background and e = 0.5 as a function 
of the spatial grid variable n. The left panel shows the inter-site solutions, while the right one shows the on-site solutions. 
The solutions to the cubic model are marked by stars, connected (as a guide to the eye) by a dashed line, while those to the 
saturable model are shown by circles, connected by a solid line. 

by means of the Newton-Raphson root-finding algorithm, as fixed points of a numerical iteration scheme on the lattice 
grid. The stars that are connected (as a guide to the eye) by dashed lines pertain to the cubic model, as opposed to 
the circles (connected by a solid line) for the saturable model. It can easily be observed that the former kinks are 
narrower and steeper, which can be attributed to the "stronger" nature of the cubic nonlincarity (in comparison with 
the weaker saturable one, which returns to a linear behavior for large amplitudes). 

We now turn to a quantitative comparison of the stability results for the on-site configurations in Fig. [21 and the 
inter-site configurations in Fig. [3] In the case of the on-site configurations, we can clearly observe the eigenvalue 
starting from i(02 ~ A2) in the saturable case and from iA± in the cubic one; this eigenvalue moves linearly (in a 
decreasing way) along the imaginary axis as e is increased, as is predicted both by the Gerschgorin estimate, and also 
remarkably accurately (see the dashed lines in the left panels of Fig. |2J) by the anti-symmetric approximation of Eqs. 
(|48|l and 1251) respectively. This eigenvalue collides, in our numerical results, with the band edge of the continuous 
spectrum for e 0.055 in the saturable model and e ss 0.077 in the cubic one; this is in remarkable agreement with the 
theoretical predictions of 0.0538 for the saturable case and 0.07735 for the cubic one, respectively given by Eqs. if^l 
and l|26|) . Notice that the latter case was the one examined previously in (2?|]. After the collision, the Hamiltonian 
Hopf bifurcation ensues, resulting into the emergence of a complex quartet of eigenvalues in the spectral plane. This 
eigenvalue approaches the spectral plane origin at A = 0, as the continuum limit of e — * 00 is approached. Here, 
it is worth briefly iterating on one of the important points of p7j for reasons of completeness. Our computations, 
performed with lattices of N = 250 lattice sites, seem to illustrate the presence of restabilization windows where the 
eigenvalue "sneaks into" the imaginary axis as e grows. However, this is a consequence of the finite computational size 
of the lattice, which leads to quantization of the relevant wavenumbers k and, as a result, to the presence of gaps in 
the continuous spectral band of Eqs. Ij44(l and 1)21(1. This is also illustrated in the top right panel of Fig. ^featuring 
(by a dashed line) a computation with a lattice four times larger than the original one, where many of the original 
gaps have been eliminated. Notice, however, that this issue (of finite size stabilization) may be one that is relevant to 
experimental situations as, e.g., in the experiments of |25| where propagation over 250 channels was reported (and, 
in fact, the generated beam was only about 25 channels wide). 

For the inter-site dark solitons, we show the relevant stability results in Fig. |3| both for the saturable (left panel) 
and the cubic (right panel) nonlinearity. The theoretical predictions of Eqs. I|46l) and (|23() are also shown by dashed 
lines and provide a fair approximation of the relevant eigenvalue, especially for small e (the perturbative result is 
not expected to be accurate for e > 0.2). An important observation here is the fact that the real eigenvalue pair, 
predicted to bifurcate immediately for e > for the inter-site configuration, remains always real up to the continuum 
limit, asymptotically approaching the spectral plane origin as e — * 00. Interestingly, this is so both for the saturable 
and for the cubic nonlinearity cases and this feature is unaffected by the lattice size (we repeated the computation of 
the saturable case with a four times larger lattice, only to find an identical pair of real eigenvalues). 

The above result also highlights a noteworthy difference between the saturable model in its defocusing and focusing 
[2(| forms. The focusing version of the model has the remarkable feature of stability alternation between on-site and 
inter-site configurations (vanishing of the so-called Peierls-Nabarro barrier) as one of its most important characteristics. 
This is apparently occurring 35] through an exchange of a real eigenvalue pair (pertaining to one of the configurations) 
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FIG. 2: (Color online) The figure shows stability results for the on-site dark soliton. The top panels pertain to the saturable 
model, while the bottom ones to the cubic model. The left panels show the trajectory of the imaginary eigenvalue as a function 
of e, starting at i(/?2 — A2) and at iAi in the two models (the maximal imaginary eigenvalue) until the point of collision e cr 
with the upper band edge of the continuous spectrum. The dashed lines in both cases show the theoretical predictions of Eqs. 
148> and 12511 respectively. The right panels show the real part of the relevant eigenvalue, which is zero before -and non-zero 
after- the relevant collision due to the ensuing Hamiltonian Hopf bifurcation leading to the exit of an eigenvalue quartet. In 
the top (saturable) case, notice in addition to the solid line (computation on a lattice with 250 sites) the dashed line which 
denotes the same computation on a lattice with 1000 sites. 




FIG. 3: (Color online) The figure shows the stability results for the inter-site dark soliton. The left and right panels correspond 
to the saturable and cubic nonlinearity cases, respectively. The solid lines refer to the numerical results, while the dashed ones 
to the theoretical predictions of Eqs. 1461 and (12311 . respectively. Notice that the configurations are unstable Ve in both cases, 
due to the relevant real eigenvalue pair. 
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FIG. 4: The figure shows a semi-log plot of the difference AG = G B — G A of the grand canonical energies G A and G B [see 
Eq. 1511 1 of the on-site and inter-site dark solitons, respectively, as a function of e. The left and right panels correspond to the 
cases of the cubic and the saturable nonlinearity, respectively. 



and of an imaginary pair (pertaining to the other) at "exceptional" values of e, marking these so-called transparency 
windows. Interestingly, this feature appears to be absent from the defocusing version of the model, possibly because 
of the structurally different nature of the spectrum and the presence of the Hamiltonian Hopf bifurcation, which 
does not allow such an exchange of eigenvalue pairs to occur at any finite value of e. We corroborate this point by 
considering the grand canonical free energy defined as 



(51) 



where j = 1,2 for the cubic and saturable nonlinearity cases, respectively, N — \u n \ is the power (or the number 
of atoms in the BEC context) and Hj denotes the Hamiltonian of the respective models, which is given by 



( \2 I 0t\ 

(U n+ 1-U n ) + —\U r , 



(52) 



for the cubic nonlinearity, and 



H 2 = 



[(un+i - u n f - /3 2 ln(l + |u„| 2 )] 



(53) 



for the saturable nonlinearity. G (and, in particular, its differences between on-site and inter-site configurations) was 
shown in 36] to be the proper indicator of the location of the exchange of stability /transparency windows in the 
focusing version of the model. Therefore, in Fig. 0|we show the difference AG = G B — G A of the grand canonical 
free energies G A and G B of the on-site and inter-site dark solitons, as a function of e in the two models. This 
quantity has the sense of a "pinning energy barrier" and is a measure of the well-known Peirls-Nabarro barrier. As 
is clearly seen in Fig. 0] AG is a monotonically decreasing function of e in both the cubic and saturable nonlinearity 
cases (shown, respectively, in the left and right panels of Fig. HJ; note that a similar result has been obtained in 
29] for the cubic case. This result is a clear indication that, in contrast to the saturable self-focusing nonlinearity 
case, in the self-defocusing regime the saturable lattice possesses no transparency points. It remains as an interesting 
open question whether such a stability alternation could be achievable in other defocusing lattice models (such as, 
e.g., possibly a cubic-quintic equation), especially because these transparency windows have been associated with 
important travelling solution properties of these lattices |36) . 

We also briefly compare the growth rates of the instabilities in the different models. In the saturable case and for 
the on-site configuration, the maximal growth rate (i.e., the maximal real part of the most unstable eigenvalue) is 
~ 0.0339 and occurs for e w 0.168. On the other hand, for the cubic case it is w 0.1199 (for e w 0.322). For the 
inter-site case, the maximal rate is 0.1972 for the saturable case (and occurs for e « 0.138), while it is 0.5163 for the 
cubic one (for e w 0.234). This illustrates a general feature that we have observed (and which is justified on the basis of 
the nature of the nonlinearity) , namely that the more unstable situations emerge for the "stronger" cubic nonlinearity 
in comparison with its saturable counterpart. This is also evinced by the dynamical evolution of the instability for 
a given e(= 0.5 in the case of Fig. [SJ. In particular, the figure shows the spatio-temporal contour plot evolution of 
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FIG. 5: (Color online) The evolution of the instability for the saturable on-site (top left), saturable inter-site (top right), cubic 
on-site (bottom left) and cubic inter-site (bottom right) dark solitons. The dynamics is shown by a space-time contour plot of 
the square modulus |w n (t)| 2 of the field. 

the square modulus |u n | 2 (£) of the field. It is clear that the dark soliton of the saturable model (top panels) with 
the weaker growth rates becomes unstable for longer times than its cubic counterpart. Moreover, the inter-site dark 
soliton (with the real eigenvalue and the -generally- larger growth rates than the eigenvalue quartet of the on-site 
case) of the right panels becomes unstable more rapidly than the on-site soliton of the left panels. The instability 
growth rates for this particular value of e are 0.0051 for the on-site saturable case (top left), 0.0939 for the inter-site 
saturable (top right), 0.0774 for the on-site cubic (bottom left) and 0.3944 for the inter-site cubic (bottom right). The 
ordering of the growth rates can be readily (and inversely) associated with the time it takes for the corresponding 
structures to become unstable. Finally, we note in passing that the oscillatory instability is manifested through an 
oscillation of the center, as is expected based on the relevant complex eigenvalue and its associated oscillatory growth. 
At the same time, the real eigenvalue pairs lead to a uni-directional propagation of the dark soliton as a result of the 
instability. This corroborates the results of [2!j for the cubic case. 

IV. CONCLUSIONS AND FUTURE CHALLENGES 

In the present paper, motivated by recent experiments in waveguide arrays featuring a saturable self-defocusing 
nonlinearity |25j . as well as earlier experimental works |24j in arrays exhibiting cubic (Kerr) nonlinearities, we have 
examined and compared two correspondingly relevant nonlinear dynamical lattice models, namely the one with a 
saturable nonlinearity, as well as the intensively examined one with the cubic nonlinearity. The objective of the study 
was to obtain an analytical, to the extent possible, understanding of the features of dark solitons, reported in earlier 
works, either computational or experimental ones (such as those cited above). 

Earlier reported features included, among others, the stable nature of on-site dark solitons versus the unstable nature 
of inter-site ones. We rationalized this on the basis of an immediately unstable (off of the anti-continuum limit) real 
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eigenvalue pair in the latter case, versus an eigenvalue on the imaginary axis, per the zero site, in the former case. 
However, it was illustrated (corroborating and expanding the earlier results of |27ll29| '). how this imaginary eigenvalue 
leads to a Hamiltonian-Hopf bifurcation eventually unveiling an oscillatorily unstable behavior in the dynamics of the 
corresponding on-site dark solitons. 

These earlier, including very recent, results also indicated (see e.g. |25|) an apparent stabilization of even inter-site 
modes in the saturable case (in comparison with the cubic one). We illustrated that while this does not appear to be 
due to a true stabilization of the solutions, it is justified by the much weaker instability growth rates of the saturable 
model, compared to its more strongly nonlinear cubic sibling. 

In addition to the above, we highlighted the usefulness of the anti-continuum limit and of perturbative expansions 
off of it, in order to obtain not only qualitative but also quantitative insight regarding the relevant structures and 
their form, as well as, more importantly, their stability. 

These techniques can be readily applied to other defocusing lattices as well and it would be of interest to do so. 
This is especially so in light of the fact that the present models did not feature an apparent exchange of stability such 
as the one showcased by the focusing analog of the saturable model. It would be interesting to examine whether other 
models, such as ones with competing nonlinearities in the form of, e.g., the cubic-quintic model might present such 
stability alternations (possibly using the extra "knob" of manipulating the strength of the quintic term). Of equal 
interest would be the expansion of considerations developed herein into higher-dimensional defocusing lattices and 
to a potential identification of the natural multi-dimensional generalization of the dark soliton solutions presented 
herein in the form of a vortex type structure. Such works are currently in progress and will be reported in future 
publications. 

Acknowledgements. PGK gratefully acknowledges support from NSF-DMS-0204585, NSF-DMS-0505663 and 
NSF-CAREER. 



APPENDIX A: PERTURBATION THEORY USING FORMAL EXPANSION 



In this appendix, as advertised above, we present another method of deriving the approximation to the key eigen- 
values of discrete dark solitons both in cubic and saturable nonlinearities and for both inter-site and on-site cases. 
Here we will use the mathematically familiar approach of formal perturbation expansions. 



1. Cubic nonlinear ity 

First, let us consider the simpler case, i.e., the cubic nonlinearity case. 

To perform the linear stability analysis around the discrete solitary waves, we introduce the following linearization 
ansatz 

u n = v n + SC„. 

Substituting into ||TJ yields the following linearized equation to Q{$) 

iC n = -eA 2 C n + fa (2\v n \ 2 C n + v 2 n C*) - A x C n . (Al) 
Decomposing C n (t) =r) n + i£ n and assuming that v n is real, eq. (|A1() gives (see, e.g., 

£)-(-/Vo w ) (£)-*(£)■ (A2 » 

where the operators £-(e) and £+(e) are defined as £-(e) = — eA 2 + 0iv\ — Ai and £+(e) = — eA 2 + 3/3i«^ — Ai. 
The stability of v n is then determined by the eigenvalues of Ti. 

Let the eigenvalues of Ti be denoted by A, which implies that v n is stable if A r = 0. Because (|A2() is linear, we can 
eliminate one of the eigenvectors, for instance £„, from which we obtain the following eigenvalue problem 

£_(e)£ + (e)?7 n = -A 2 r/„ = Zrj n . (A3) 

As before, we expand the eigenvector rj n and the eigenvalue 5 as 

Vn = r,P + + 0{e% E = S<°> + eZ^ + {e 2 ). 
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Substituting into Eq. (|A3|) and identifying coefficients for consecutive powers of the small parameter e yields 



£_(0)£ + (0)-HW> 
£_(0)£+(0)-S(°) 



W = o, 

€ ] = f, 



with 



/ 



(A 2 - 2(3 lV ^v^)£ + (0) + £_(0)(A a - 6«4%«) + 5 



(A4) 
(A5) 

(A6) 



First, let us consider the order 0(1) equation (|A4|I . One can do a simple analysis as above to show that there is one 
eigenvalue, i.e., — for the inter-site configuration and two eigenvalues = and S^ ^ = Af for the on-site 
one. The zero eigenvalue has infinite multiplicity and is related to the continuous spectrum, as discussed previously. 

For the inter-site configuration, there is an eigenvalue bifurcating from the continuous spectrum as soon as the 
coupling is turned on. Therefore, this zero eigenvalue is the crucial eigenvalue for its stability. The normalized 



eigenvector of this eigenvalue is rj n =1/ a/2, for 



0, 1 and r/n 



otherwise. For the on-site configuration, the 



A^ with the normalized eigenvector r) n °^ = 1, for n 



and r]n^ 







crucial eigenvalue for the stability is 
otherwise. 

The continuation of the critical eigenvalue when the coupling e is turned on can be calculated from Eq. (|A5() . Due to 
the fact that the corresponding eigenvector is zero almost everywhere, we only need to consider the site with non-zero 
component eigenvector, i.e. n = 0, 1 for the inter-site and n = for the on-site. 

It is simple to show that the solvability condition of Eq. I|A5(1 using, e.g., the Fredholm alternative requires / = 
from which one immeadiately obtains that = — 2Ai for the inter-site and E^ 1 ) = — 4Ai for the inter-site. 

Hence, the critical eigenvalue is: 



X = ±^/2A 1 e + 0(e 2 ) 



for the inter-site configuration and 



A = ±i 



A? 



4Ae + C(e 2 ) 



(A7) 



(A8) 



for the on-site one. 

These two results agree with the leading order calculations of section II. A. In particular, Eq. (|A7(I can be immedi- 
ately seen to agree with the leading order prediction of Eq. I|23|l . On the other hand, as regards Eq. (|A8|) . a Taylor 
expansion to leading order yields A = ±i(Ai — 2e), again in agreement with the findings of Eq. Q25J1. 



2. Saturable nonlinearity 



For the saturable nonlinearity, substituting the same spectral ansatz u n — v n + SC n into l|27() yields the following 
linearized equation to 0(5) 



iC n = -eA 2 C n 



1 + |v„p 



VnjVnC* + V*C n ) 

i + KI 2 



-A 2 C„. 



(A9) 



With this nonlinearity, the operator £_ (e) and £+(e) are now defined as (e) = — eA 2 + ^(l + v^) /(1 + |u„| 2 ) 2 — A9 
and C+(e) = -eA 2 + 2 (1 - v 2 n )/{\ + \v n \ 2 ) 2 ~ A 2 . 



Again, by using S — S^ ^ + e^/ 1 -* and rj n — ?y„ + erf n , we obtain the same equations as Eqs. (|A4IA5p with the 
inhomogeneity term / now given by 



(A 2 - k ln )C+(0) + £_(0)(A 2 - k 2n ) + S« 



„(0) 
'In ' 



(A10) 



with 



kin — 



2(3 2 v^v^ 



k 2n = - 



2/3 2 ^ 0) ^ 1) 



(l + og ,) ) 



d+i4 0) ) 



2(1~,1 0) ) 
(l + ^° )2 ) 
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The stability of the static solutions of Eqs. (|31|l - (|35|l and l|39|) - <|42[) can then be determined using exactly the same 
procedures as above. One can simply calculate that the critical eigenvalue is 



for the on-site configuration. 

The eigenvalue prediction from Eq. (|A11|I can be immediately seen to be identical to the leading order prediction 
of Eq. (SB) . On the other hand, Taylor expanding the result of Eq. i|A12|) . we obtain A = ±i(/3 2 — A 2 — 2e), which 
once again is identical to the prediction of Eq. (|48|l . 
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